<!-- HTML header for doxygen 1.8.6-->
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<meta name="generator" content="Doxygen 1.8.13"/>
<title>OpenCV: Reading Geospatial Raster files with GDAL</title>
<link href="../../opencv.ico" rel="shortcut icon" type="image/x-icon" />
<link href="../../tabs.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="../../jquery.js"></script>
<script type="text/javascript" src="../../dynsections.js"></script>
<script type="text/javascript" src="../../tutorial-utils.js"></script>
<link href="../../search/search.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="../../search/searchdata.js"></script>
<script type="text/javascript" src="../../search/search.js"></script>
<script type="text/x-mathjax-config">
  MathJax.Hub.Config({
    extensions: ["tex2jax.js", "TeX/AMSmath.js", "TeX/AMSsymbols.js"],
    jax: ["input/TeX","output/HTML-CSS"],
});
//<![CDATA[
MathJax.Hub.Config(
{
  TeX: {
      Macros: {
          matTT: [ "\\[ \\left|\\begin{array}{ccc} #1 & #2 & #3\\\\ #4 & #5 & #6\\\\ #7 & #8 & #9 \\end{array}\\right| \\]", 9],
          fork: ["\\left\\{ \\begin{array}{l l} #1 & \\mbox{#2}\\\\ #3 & \\mbox{#4}\\\\ \\end{array} \\right.", 4],
          forkthree: ["\\left\\{ \\begin{array}{l l} #1 & \\mbox{#2}\\\\ #3 & \\mbox{#4}\\\\ #5 & \\mbox{#6}\\\\ \\end{array} \\right.", 6],
          forkfour: ["\\left\\{ \\begin{array}{l l} #1 & \\mbox{#2}\\\\ #3 & \\mbox{#4}\\\\ #5 & \\mbox{#6}\\\\ #7 & \\mbox{#8}\\\\ \\end{array} \\right.", 8],
          vecthree: ["\\begin{bmatrix} #1\\\\ #2\\\\ #3 \\end{bmatrix}", 3],
          vecthreethree: ["\\begin{bmatrix} #1 & #2 & #3\\\\ #4 & #5 & #6\\\\ #7 & #8 & #9 \\end{bmatrix}", 9],
          cameramatrix: ["#1 = \\begin{bmatrix} f_x & 0 & c_x\\\\ 0 & f_y & c_y\\\\ 0 & 0 & 1 \\end{bmatrix}", 1],
          distcoeffs: ["(k_1, k_2, p_1, p_2[, k_3[, k_4, k_5, k_6 [, s_1, s_2, s_3, s_4[, \\tau_x, \\tau_y]]]]) \\text{ of 4, 5, 8, 12 or 14 elements}"],
          distcoeffsfisheye: ["(k_1, k_2, k_3, k_4)"],
          hdotsfor: ["\\dots", 1],
          mathbbm: ["\\mathbb{#1}", 1],
          bordermatrix: ["\\matrix{#1}", 1]
      }
  }
}
);
//]]>
</script><script type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/MathJax.js"></script>
<link href="../../doxygen.css" rel="stylesheet" type="text/css" />
<link href="../../stylesheet.css" rel="stylesheet" type="text/css"/>
</head>
<body>
<div id="top"><!-- do not remove this div, it is closed by doxygen! -->
<div id="titlearea">
<!--#include virtual="/google-search.html"-->
<table cellspacing="0" cellpadding="0">
 <tbody>
 <tr style="height: 56px;">
  <td id="projectlogo"><img alt="Logo" src="../../opencv-logo-small.png"/></td>
  <td style="padding-left: 0.5em;">
   <div id="projectname">OpenCV
   &#160;<span id="projectnumber">4.5.2</span>
   </div>
   <div id="projectbrief">Open Source Computer Vision</div>
  </td>
 </tr>
 </tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.8.13 -->
<script type="text/javascript">
var searchBox = new SearchBox("searchBox", "../../search",false,'Search');
</script>
<script type="text/javascript" src="../../menudata.js"></script>
<script type="text/javascript" src="../../menu.js"></script>
<script type="text/javascript">
$(function() {
  initMenu('../../',true,false,'search.php','Search');
  $(document).ready(function() { init_search(); });
});
</script>
<div id="main-nav"></div>
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
     onmouseover="return searchBox.OnSearchSelectShow()"
     onmouseout="return searchBox.OnSearchSelectHide()"
     onkeydown="return searchBox.OnSearchSelectKey(event)">
</div>

<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0" 
        name="MSearchResults" id="MSearchResults">
</iframe>
</div>

<div id="nav-path" class="navpath">
  <ul>
<li class="navelem"><a class="el" href="../../d9/df8/tutorial_root.html">OpenCV Tutorials</a></li><li class="navelem"><a class="el" href="../../de/d3d/tutorial_table_of_content_app.html">Application utils (highgui, imgcodecs, videoio modules)</a></li>  </ul>
</div>
</div><!-- top -->
<div class="header">
  <div class="headertitle">
<div class="title">Reading Geospatial Raster files with GDAL </div>  </div>
</div><!--header-->
<div class="contents">
<div class="textblock"><p><b>Prev Tutorial:</b> <a class="el" href="../../da/d6a/tutorial_trackbar.html">Adding a Trackbar to our applications!</a></p>
<p><b>Next Tutorial:</b> <a class="el" href="../../d5/dc4/tutorial_video_input_psnr_ssim.html">Video Input with OpenCV and similarity measurement</a></p>
<table class="doxtable">
<tr>
<th align="right"></th><th align="left"></th></tr>
<tr>
<td align="right">Original author </td><td align="left">Marvin Smith </td></tr>
<tr>
<td align="right">Compatibility </td><td align="left">OpenCV &gt;= 3.0 </td></tr>
</table>
<p>Geospatial raster data is a heavily used product in Geographic Information Systems and Photogrammetry. Raster data typically can represent imagery and Digital Elevation Models (DEM). The standard library for loading GIS imagery is the Geographic Data Abstraction Library <a href="http://www.gdal.org">(GDAL)</a>. In this example, we will show techniques for loading GIS raster formats using native OpenCV functions. In addition, we will show some an example of how OpenCV can use this data for novel and interesting purposes.</p>
<h2>Goals </h2>
<p>The primary objectives for this tutorial:</p>
<ul>
<li>How to use OpenCV <a class="el" href="../../d4/da8/group__imgcodecs.html#imread">imread</a> to load satellite imagery.</li>
<li>How to use OpenCV <a class="el" href="../../d4/da8/group__imgcodecs.html#imread">imread</a> to load SRTM Digital Elevation Models</li>
<li>Given the corner coordinates of both the image and DEM, correlate the elevation data to the image to find elevations for each pixel.</li>
<li>Show a basic, easy-to-implement example of a terrain heat map.</li>
<li>Show a basic use of DEM data coupled with ortho-rectified imagery.</li>
</ul>
<p>To implement these goals, the following code takes a Digital Elevation Model as well as a GeoTiff image of San Francisco as input. The image and DEM data is processed and generates a terrain heat map of the image as well as labels areas of the city which would be affected should the water level of the bay rise 10, 50, and 100 meters.</p>
<h2>Code </h2>
<div class="fragment"><div class="line"><span class="comment">/*</span></div><div class="line"><span class="comment"> * gdal_image.cpp -- Load GIS data into OpenCV Containers using the Geospatial Data Abstraction Library</span></div><div class="line"><span class="comment">*/</span></div><div class="line"></div><div class="line"><span class="comment">// OpenCV Headers</span></div><div class="line"><span class="preprocessor">#include &quot;<a class="code" href="../../d0/d9c/core_2include_2opencv2_2core_8hpp.html">opencv2/core.hpp</a>&quot;</span></div><div class="line"><span class="preprocessor">#include &quot;<a class="code" href="../../d1/d4f/imgproc_2include_2opencv2_2imgproc_8hpp.html">opencv2/imgproc.hpp</a>&quot;</span></div><div class="line"><span class="preprocessor">#include &quot;<a class="code" href="../../d4/dd5/highgui_8hpp.html">opencv2/highgui.hpp</a>&quot;</span></div><div class="line"></div><div class="line"><span class="comment">// C++ Standard Libraries</span></div><div class="line"><span class="preprocessor">#include &lt;cmath&gt;</span></div><div class="line"><span class="preprocessor">#include &lt;iostream&gt;</span></div><div class="line"><span class="preprocessor">#include &lt;stdexcept&gt;</span></div><div class="line"><span class="preprocessor">#include &lt;vector&gt;</span></div><div class="line"></div><div class="line"><span class="keyword">using namespace </span><a class="code" href="../../d8/dcc/namespacestd.html">std</a>;</div><div class="line"></div><div class="line"><span class="comment">// define the corner points</span></div><div class="line"><span class="comment">//    Note that GDAL library can natively determine this</span></div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> tl( -122.441017, 37.815664 );</div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> tr( -122.370919, 37.815311 );</div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> bl( -122.441533, 37.747167 );</div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> br( -122.3715,   37.746814 );</div><div class="line"></div><div class="line"><span class="comment">// determine dem corners</span></div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> dem_bl( -122.0, 38);</div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> dem_tr( -123.0, 37);</div><div class="line"></div><div class="line"><span class="comment">// range of the heat map colors</span></div><div class="line">std::vector&lt;std::pair&lt;cv::Vec3b,double&gt; &gt; color_range;</div><div class="line"></div><div class="line"></div><div class="line"><span class="comment">// List of all function prototypes</span></div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> lerp( <span class="keyword">const</span> <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a>&amp;, <span class="keyword">const</span> <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a>&amp;, <span class="keyword">const</span> <span class="keywordtype">double</span>&amp; );</div><div class="line"></div><div class="line"><a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a> get_dem_color( <span class="keyword">const</span> <span class="keywordtype">double</span>&amp; );</div><div class="line"></div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> world2dem( <span class="keyword">const</span> <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a>&amp;, <span class="keyword">const</span> <a class="code" href="../../d6/d50/classcv_1_1Size__.html">cv::Size</a>&amp;);</div><div class="line"></div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> pixel2world( <span class="keyword">const</span> <span class="keywordtype">int</span>&amp;, <span class="keyword">const</span> <span class="keywordtype">int</span>&amp;, <span class="keyword">const</span> <a class="code" href="../../d6/d50/classcv_1_1Size__.html">cv::Size</a>&amp; );</div><div class="line"></div><div class="line"><span class="keywordtype">void</span> add_color( <a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>&amp; pix, <span class="keyword">const</span> <a class="code" href="../../d1/d1b/group__core__hal__interface.html#ga65f85814a8290f9797005d3b28e7e5fc">uchar</a>&amp; b, <span class="keyword">const</span> <a class="code" href="../../d1/d1b/group__core__hal__interface.html#ga65f85814a8290f9797005d3b28e7e5fc">uchar</a>&amp; g, <span class="keyword">const</span> <a class="code" href="../../d1/d1b/group__core__hal__interface.html#ga65f85814a8290f9797005d3b28e7e5fc">uchar</a>&amp; r );</div><div class="line"></div><div class="line"></div><div class="line"></div><div class="line"><span class="comment">/*</span></div><div class="line"><span class="comment"> * Linear Interpolation</span></div><div class="line"><span class="comment"> * p1 - Point 1</span></div><div class="line"><span class="comment"> * p2 - Point 2</span></div><div class="line"><span class="comment"> * t  - Ratio from Point 1 to Point 2</span></div><div class="line"><span class="comment">*/</span></div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> lerp( <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> <span class="keyword">const</span>&amp; p1, <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> <span class="keyword">const</span>&amp; p2, <span class="keyword">const</span> <span class="keywordtype">double</span>&amp; t ){</div><div class="line">    <span class="keywordflow">return</span> <a class="code" href="../../dc/d84/group__core__basic.html#ga32986d8c6fd05949688352abd112509f">cv::Point2d</a>( ((1-t)*p1.<a class="code" href="../../db/d4e/classcv_1_1Point__.html#a4c96fa7bdbfe390be5ed356edb274ff3">x</a>) + (t*p2.<a class="code" href="../../db/d4e/classcv_1_1Point__.html#a4c96fa7bdbfe390be5ed356edb274ff3">x</a>),</div><div class="line">                        ((1-t)*p1.<a class="code" href="../../db/d4e/classcv_1_1Point__.html#a157337197338ff199e5df1a393022f15">y</a>) + (t*p2.<a class="code" href="../../db/d4e/classcv_1_1Point__.html#a157337197338ff199e5df1a393022f15">y</a>));</div><div class="line">}</div><div class="line"></div><div class="line"><span class="comment">/*</span></div><div class="line"><span class="comment"> * Interpolate Colors</span></div><div class="line"><span class="comment">*/</span></div><div class="line"><span class="keyword">template</span> &lt;<span class="keyword">typename</span> DATATYPE, <span class="keywordtype">int</span> N&gt;</div><div class="line"><a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec&lt;DATATYPE,N&gt;</a> lerp( <a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec&lt;DATATYPE,N&gt;</a> <span class="keyword">const</span>&amp; minColor,</div><div class="line">                          <a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec&lt;DATATYPE,N&gt;</a> <span class="keyword">const</span>&amp; maxColor,</div><div class="line">                          <span class="keywordtype">double</span> <span class="keyword">const</span>&amp; t ){</div><div class="line"></div><div class="line">    <a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec&lt;DATATYPE,N&gt;</a> output;</div><div class="line">    <span class="keywordflow">for</span>( <span class="keywordtype">int</span> i=0; i&lt;N; i++ ){</div><div class="line">        output[i] = (<a class="code" href="../../d1/d1b/group__core__hal__interface.html#ga65f85814a8290f9797005d3b28e7e5fc">uchar</a>)(((1-t)*minColor[i]) + (t * maxColor[i]));</div><div class="line">    }</div><div class="line">    <span class="keywordflow">return</span> output;</div><div class="line">}</div><div class="line"></div><div class="line"><span class="comment">/*</span></div><div class="line"><span class="comment"> * Compute the dem color</span></div><div class="line"><span class="comment">*/</span></div><div class="line"><a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a> get_dem_color( <span class="keyword">const</span> <span class="keywordtype">double</span>&amp; elevation ){</div><div class="line"></div><div class="line">    <span class="comment">// if the elevation is below the minimum, return the minimum</span></div><div class="line">    <span class="keywordflow">if</span>( elevation &lt; color_range[0].second ){</div><div class="line">        <span class="keywordflow">return</span> color_range[0].first;</div><div class="line">    }</div><div class="line">    <span class="comment">// if the elevation is above the maximum, return the maximum</span></div><div class="line">    <span class="keywordflow">if</span>( elevation &gt; color_range.back().second ){</div><div class="line">        <span class="keywordflow">return</span> color_range.back().first;</div><div class="line">    }</div><div class="line"></div><div class="line">    <span class="comment">// otherwise, find the proper starting index</span></div><div class="line">    <span class="keywordtype">int</span> idx=0;</div><div class="line">    <span class="keywordtype">double</span> t = 0;</div><div class="line">    <span class="keywordflow">for</span>( <span class="keywordtype">int</span> x=0; x&lt;(int)(color_range.size()-1); x++ ){</div><div class="line"></div><div class="line">        <span class="comment">// if the current elevation is below the next item, then use the current</span></div><div class="line">        <span class="comment">// two colors as our range</span></div><div class="line">        <span class="keywordflow">if</span>( elevation &lt; color_range[x+1].second ){</div><div class="line">            idx=x;</div><div class="line">            t = (color_range[x+1].second - elevation)/</div><div class="line">                (color_range[x+1].second - color_range[x].second);</div><div class="line"></div><div class="line">            <span class="keywordflow">break</span>;</div><div class="line">        }</div><div class="line">    }</div><div class="line"></div><div class="line">    <span class="comment">// interpolate the color</span></div><div class="line">    <span class="keywordflow">return</span> lerp( color_range[idx].first, color_range[idx+1].first, t);</div><div class="line">}</div><div class="line"></div><div class="line"><span class="comment">/*</span></div><div class="line"><span class="comment"> * Given a pixel coordinate and the size of the input image, compute the pixel location</span></div><div class="line"><span class="comment"> * on the DEM image.</span></div><div class="line"><span class="comment">*/</span></div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> world2dem( <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> <span class="keyword">const</span>&amp; coordinate, <span class="keyword">const</span> <a class="code" href="../../d6/d50/classcv_1_1Size__.html">cv::Size</a>&amp; dem_size   ){</div><div class="line"></div><div class="line"></div><div class="line">    <span class="comment">// relate this to the dem points</span></div><div class="line">    <span class="comment">// ASSUMING THAT DEM DATA IS ORTHORECTIFIED</span></div><div class="line">    <span class="keywordtype">double</span> demRatioX = ((dem_tr.x - coordinate.<a class="code" href="../../db/d4e/classcv_1_1Point__.html#a4c96fa7bdbfe390be5ed356edb274ff3">x</a>)/(dem_tr.x - dem_bl.x));</div><div class="line">    <span class="keywordtype">double</span> demRatioY = 1-((dem_tr.y - coordinate.<a class="code" href="../../db/d4e/classcv_1_1Point__.html#a157337197338ff199e5df1a393022f15">y</a>)/(dem_tr.y - dem_bl.y));</div><div class="line"></div><div class="line">    <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> output;</div><div class="line">    output.<a class="code" href="../../db/d4e/classcv_1_1Point__.html#a4c96fa7bdbfe390be5ed356edb274ff3">x</a> = demRatioX * dem_size.<a class="code" href="../../d6/d50/classcv_1_1Size__.html#abfe0367b32c407ddccf5ddf92667c73d">width</a>;</div><div class="line">    output.<a class="code" href="../../db/d4e/classcv_1_1Point__.html#a157337197338ff199e5df1a393022f15">y</a> = demRatioY * dem_size.<a class="code" href="../../d6/d50/classcv_1_1Size__.html#a1d289dce6b5d8006a54f3ee0259fc545">height</a>;</div><div class="line"></div><div class="line">    <span class="keywordflow">return</span> output;</div><div class="line">}</div><div class="line"></div><div class="line"><span class="comment">/*</span></div><div class="line"><span class="comment"> * Convert a pixel coordinate to world coordinates</span></div><div class="line"><span class="comment">*/</span></div><div class="line"><a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> pixel2world( <span class="keyword">const</span> <span class="keywordtype">int</span>&amp; x, <span class="keyword">const</span> <span class="keywordtype">int</span>&amp; y, <span class="keyword">const</span> <a class="code" href="../../d6/d50/classcv_1_1Size__.html">cv::Size</a>&amp; <a class="code" href="../../df/d5b/namespacecv_1_1gapi_1_1streaming.html#a0a915e69f4cc8284293e40fc9ffbf157">size</a> ){</div><div class="line"></div><div class="line">    <span class="comment">// compute the ratio of the pixel location to its dimension</span></div><div class="line">    <span class="keywordtype">double</span> rx = (double)x / size.<a class="code" href="../../d6/d50/classcv_1_1Size__.html#abfe0367b32c407ddccf5ddf92667c73d">width</a>;</div><div class="line">    <span class="keywordtype">double</span> ry = (<span class="keywordtype">double</span>)y / size.<a class="code" href="../../d6/d50/classcv_1_1Size__.html#a1d289dce6b5d8006a54f3ee0259fc545">height</a>;</div><div class="line"></div><div class="line">    <span class="comment">// compute LERP of each coordinate</span></div><div class="line">    <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> rightSide = lerp(tr, br, ry);</div><div class="line">    <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> leftSide  = lerp(tl, bl, ry);</div><div class="line"></div><div class="line">    <span class="comment">// compute the actual Lat/Lon coordinate of the interpolated coordinate</span></div><div class="line">    <span class="keywordflow">return</span> lerp( leftSide, rightSide, rx );</div><div class="line">}</div><div class="line"></div><div class="line"><span class="comment">/*</span></div><div class="line"><span class="comment"> * Add color to a specific pixel color value</span></div><div class="line"><span class="comment">*/</span></div><div class="line"><span class="keywordtype">void</span> add_color( <a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>&amp; pix, <span class="keyword">const</span> <a class="code" href="../../d1/d1b/group__core__hal__interface.html#ga65f85814a8290f9797005d3b28e7e5fc">uchar</a>&amp; b, <span class="keyword">const</span> <a class="code" href="../../d1/d1b/group__core__hal__interface.html#ga65f85814a8290f9797005d3b28e7e5fc">uchar</a>&amp; g, <span class="keyword">const</span> <a class="code" href="../../d1/d1b/group__core__hal__interface.html#ga65f85814a8290f9797005d3b28e7e5fc">uchar</a>&amp; r ){</div><div class="line"></div><div class="line">    <span class="keywordflow">if</span>( pix[0] + b &lt; 255 &amp;&amp; pix[0] + b &gt;= 0 ){ pix[0] += b; }</div><div class="line">    <span class="keywordflow">if</span>( pix[1] + g &lt; 255 &amp;&amp; pix[1] + g &gt;= 0 ){ pix[1] += g; }</div><div class="line">    <span class="keywordflow">if</span>( pix[2] + r &lt; 255 &amp;&amp; pix[2] + r &gt;= 0 ){ pix[2] += r; }</div><div class="line">}</div><div class="line"></div><div class="line"></div><div class="line"><span class="comment">/*</span></div><div class="line"><span class="comment"> * Main Function</span></div><div class="line"><span class="comment">*/</span></div><div class="line"><span class="keywordtype">int</span> main( <span class="keywordtype">int</span> argc, <span class="keywordtype">char</span>* argv[] ){</div><div class="line"></div><div class="line">    <span class="comment">/*</span></div><div class="line"><span class="comment">     * Check input arguments</span></div><div class="line"><span class="comment">    */</span></div><div class="line">    <span class="keywordflow">if</span>( argc &lt; 3 ){</div><div class="line">        cout &lt;&lt; <span class="stringliteral">&quot;usage: &quot;</span> &lt;&lt; argv[0] &lt;&lt; <span class="stringliteral">&quot; &lt;image_name&gt; &lt;dem_model_name&gt;&quot;</span> &lt;&lt; endl;</div><div class="line">        <span class="keywordflow">return</span> -1;</div><div class="line">    }</div><div class="line"></div><div class="line">    <span class="comment">// load the image (note that we don&#39;t have the projection information.  You will</span></div><div class="line">    <span class="comment">// need to load that yourself or use the full GDAL driver.  The values are pre-defined</span></div><div class="line">    <span class="comment">// at the top of this file</span></div><div class="line"><span class="comment"></span>    <a class="code" href="../../d3/d63/classcv_1_1Mat.html">cv::Mat</a> image = <a class="code" href="../../d4/da8/group__imgcodecs.html#ga288b8b3da0892bd651fce07b3bbd3a56">cv::imread</a>(argv[1], <a class="code" href="../../d8/d6a/group__imgcodecs__flags.html#gga61d9b0126a3e57d9277ac48327799c80a1b482916f4ffe34f0b90550f267fd896">cv::IMREAD_LOAD_GDAL</a> | <a class="code" href="../../d8/d6a/group__imgcodecs__flags.html#gga61d9b0126a3e57d9277ac48327799c80af660544735200cbe942eea09232eb822">cv::IMREAD_COLOR</a> );</div><div class="line"></div><div class="line">    <span class="comment">// load the dem model</span></div><div class="line">    <a class="code" href="../../d3/d63/classcv_1_1Mat.html">cv::Mat</a> dem = <a class="code" href="../../d4/da8/group__imgcodecs.html#ga288b8b3da0892bd651fce07b3bbd3a56">cv::imread</a>(argv[2], <a class="code" href="../../d8/d6a/group__imgcodecs__flags.html#gga61d9b0126a3e57d9277ac48327799c80a1b482916f4ffe34f0b90550f267fd896">cv::IMREAD_LOAD_GDAL</a> | <a class="code" href="../../d8/d6a/group__imgcodecs__flags.html#gga61d9b0126a3e57d9277ac48327799c80a0b486c93c25e8a0b0712681bb7254c18">cv::IMREAD_ANYDEPTH</a> );</div><div class="line"></div><div class="line">    <span class="comment">// create our output products</span></div><div class="line">    <a class="code" href="../../d3/d63/classcv_1_1Mat.html">cv::Mat</a> output_dem(   image.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#a146f8e8dda07d1365a575ab83d9828d1">size</a>(), <a class="code" href="../../d1/d1b/group__core__hal__interface.html#ga88c4cd9de76f678f33928ef1e3f96047">CV_8UC3</a> );</div><div class="line">    <a class="code" href="../../d3/d63/classcv_1_1Mat.html">cv::Mat</a> output_dem_flood(   image.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#a146f8e8dda07d1365a575ab83d9828d1">size</a>(), <a class="code" href="../../d1/d1b/group__core__hal__interface.html#ga88c4cd9de76f678f33928ef1e3f96047">CV_8UC3</a> );</div><div class="line"></div><div class="line">    <span class="comment">// for sanity sake, make sure GDAL Loads it as a signed short</span></div><div class="line">    <span class="keywordflow">if</span>( dem.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#af2d2652e552d7de635988f18a84b53e5">type</a>() != <a class="code" href="../../d1/d1b/group__core__hal__interface.html#ga3515053e24377cf079fa1027d208ad52">CV_16SC1</a> ){ <span class="keywordflow">throw</span> std::runtime_error(<span class="stringliteral">&quot;DEM image type must be CV_16SC1&quot;</span>); }</div><div class="line"></div><div class="line">    <span class="comment">// define the color range to create our output DEM heat map</span></div><div class="line">    <span class="comment">//  Pair format ( Color, elevation );  Push from low to high</span></div><div class="line">    <span class="comment">//  Note:  This would be perfect for a configuration file, but is here for a working demo.</span></div><div class="line">    color_range.push_back( std::pair&lt;cv::Vec3b,double&gt;(<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>( 188, 154,  46),   -1));</div><div class="line">    color_range.push_back( std::pair&lt;cv::Vec3b,double&gt;(<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>( 110, 220, 110), 0.25));</div><div class="line">    color_range.push_back( std::pair&lt;cv::Vec3b,double&gt;(<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>( 150, 250, 230),   20));</div><div class="line">    color_range.push_back( std::pair&lt;cv::Vec3b,double&gt;(<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>( 160, 220, 200),   75));</div><div class="line">    color_range.push_back( std::pair&lt;cv::Vec3b,double&gt;(<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>( 220, 190, 170),  100));</div><div class="line">    color_range.push_back( std::pair&lt;cv::Vec3b,double&gt;(<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>( 250, 180, 140),  200));</div><div class="line"></div><div class="line">    <span class="comment">// define a minimum elevation</span></div><div class="line">    <span class="keywordtype">double</span> minElevation = -10;</div><div class="line"></div><div class="line">    <span class="comment">// iterate over each pixel in the image, computing the dem point</span></div><div class="line">    <span class="keywordflow">for</span>( <span class="keywordtype">int</span> y=0; y&lt;image.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#abed816466c45234254d25bc59c31245e">rows</a>; y++ ){</div><div class="line">    <span class="keywordflow">for</span>( <span class="keywordtype">int</span> x=0; x&lt;image.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#aa3e5a47585c9ef6a0842556739155e3e">cols</a>; x++ ){</div><div class="line"></div><div class="line">        <span class="comment">// convert the pixel coordinate to lat/lon coordinates</span></div><div class="line">        <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> coordinate = pixel2world( x, y, image.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#a146f8e8dda07d1365a575ab83d9828d1">size</a>() );</div><div class="line"></div><div class="line">        <span class="comment">// compute the dem image pixel coordinate from lat/lon</span></div><div class="line">        <a class="code" href="../../db/d4e/classcv_1_1Point__.html">cv::Point2d</a> dem_coordinate = world2dem( coordinate, dem.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#a146f8e8dda07d1365a575ab83d9828d1">size</a>() );</div><div class="line"></div><div class="line">        <span class="comment">// extract the elevation</span></div><div class="line">        <span class="keywordtype">double</span> dz;</div><div class="line">        <span class="keywordflow">if</span>( dem_coordinate.x &gt;=    0    &amp;&amp; dem_coordinate.y &gt;=    0     &amp;&amp;</div><div class="line">            dem_coordinate.x &lt; dem.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#aa3e5a47585c9ef6a0842556739155e3e">cols</a> &amp;&amp; dem_coordinate.y &lt; dem.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#abed816466c45234254d25bc59c31245e">rows</a> ){</div><div class="line">            dz = dem.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#aa5d20fc86d41d59e4d71ae93daee9726">at</a>&lt;<span class="keywordtype">short</span>&gt;(dem_coordinate);</div><div class="line">        }<span class="keywordflow">else</span>{</div><div class="line">            dz = minElevation;</div><div class="line">        }</div><div class="line"></div><div class="line">        <span class="comment">// write the pixel value to the file</span></div><div class="line">        output_dem_flood.at&lt;<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>&gt;(y,x) = image.<a class="code" href="../../d3/d63/classcv_1_1Mat.html#aa5d20fc86d41d59e4d71ae93daee9726">at</a>&lt;<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>&gt;(y,x);</div><div class="line"></div><div class="line">        <span class="comment">// compute the color for the heat map output</span></div><div class="line">        <a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a> actualColor = get_dem_color(dz);</div><div class="line">        output_dem.at&lt;<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>&gt;(y,x) = actualColor;</div><div class="line"></div><div class="line">        <span class="comment">// show effect of a 10 meter increase in ocean levels</span></div><div class="line">        <span class="keywordflow">if</span>( dz &lt; 10 ){</div><div class="line">            add_color( output_dem_flood.at&lt;<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>&gt;(y,x), 90, 0, 0 );</div><div class="line">        }</div><div class="line">        <span class="comment">// show effect of a 50 meter increase in ocean levels</span></div><div class="line">        <span class="keywordflow">else</span> <span class="keywordflow">if</span>( dz &lt; 50 ){</div><div class="line">            add_color( output_dem_flood.at&lt;<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>&gt;(y,x), 0, 90, 0 );</div><div class="line">        }</div><div class="line">        <span class="comment">// show effect of a 100 meter increase in ocean levels</span></div><div class="line">        <span class="keywordflow">else</span> <span class="keywordflow">if</span>( dz &lt; 100 ){</div><div class="line">            add_color( output_dem_flood.at&lt;<a class="code" href="../../d6/dcf/classcv_1_1Vec.html">cv::Vec3b</a>&gt;(y,x), 0, 0, 90 );</div><div class="line">        }</div><div class="line"></div><div class="line">    }}</div><div class="line"></div><div class="line">    <span class="comment">// print our heat map</span></div><div class="line">    <a class="code" href="../../d4/da8/group__imgcodecs.html#gabbc7ef1aa2edfaa87772f1202d67e0ce">cv::imwrite</a>( <span class="stringliteral">&quot;heat-map.jpg&quot;</span>   ,  output_dem );</div><div class="line"></div><div class="line">    <span class="comment">// print the flooding effect image</span></div><div class="line">    <a class="code" href="../../d4/da8/group__imgcodecs.html#gabbc7ef1aa2edfaa87772f1202d67e0ce">cv::imwrite</a>( <span class="stringliteral">&quot;flooded.jpg&quot;</span>,  output_dem_flood);</div><div class="line"></div><div class="line">    <span class="keywordflow">return</span> 0;</div><div class="line">}</div></div><!-- fragment --><h2>How to Read Raster Data using GDAL </h2>
<p>This demonstration uses the default OpenCV imread function. The primary difference is that in order to force GDAL to load the image, you must use the appropriate flag. </p><div class="fragment"><div class="line">    <a class="code" href="../../d3/d63/classcv_1_1Mat.html">cv::Mat</a> image = <a class="code" href="../../d4/da8/group__imgcodecs.html#ga288b8b3da0892bd651fce07b3bbd3a56">cv::imread</a>(argv[1], <a class="code" href="../../d8/d6a/group__imgcodecs__flags.html#gga61d9b0126a3e57d9277ac48327799c80a1b482916f4ffe34f0b90550f267fd896">cv::IMREAD_LOAD_GDAL</a> | <a class="code" href="../../d8/d6a/group__imgcodecs__flags.html#gga61d9b0126a3e57d9277ac48327799c80af660544735200cbe942eea09232eb822">cv::IMREAD_COLOR</a> );</div></div><!-- fragment --><p>When loading digital elevation models, the actual numeric value of each pixel is essential and cannot be scaled or truncated. For example, with image data a pixel represented as a double with a value of 1 has an equal appearance to a pixel which is represented as an unsigned character with a value of 255. With terrain data, the pixel value represents the elevation in meters. In order to ensure that OpenCV preserves the native value, use the GDAL flag in imread with the ANYDEPTH flag. </p><div class="fragment"><div class="line">    <span class="comment">// load the dem model</span></div><div class="line">    <a class="code" href="../../d3/d63/classcv_1_1Mat.html">cv::Mat</a> dem = <a class="code" href="../../d4/da8/group__imgcodecs.html#ga288b8b3da0892bd651fce07b3bbd3a56">cv::imread</a>(argv[2], <a class="code" href="../../d8/d6a/group__imgcodecs__flags.html#gga61d9b0126a3e57d9277ac48327799c80a1b482916f4ffe34f0b90550f267fd896">cv::IMREAD_LOAD_GDAL</a> | <a class="code" href="../../d8/d6a/group__imgcodecs__flags.html#gga61d9b0126a3e57d9277ac48327799c80a0b486c93c25e8a0b0712681bb7254c18">cv::IMREAD_ANYDEPTH</a> );</div></div><!-- fragment --><p>If you know beforehand the type of DEM model you are loading, then it may be a safe bet to test the Mat::type() or <a class="el" href="../../d0/daf/group__projection.html#ga8b5e5413af41caf3f2a0bbf4cdf517a8" title="Returns the depth of a point transformed by a rigid transform. ">Mat::depth()</a> using an assert or other mechanism. NASA or DOD specification documents can provide the input types for various elevation models. The major types, SRTM and DTED, are both signed shorts.</p>
<h2>Notes </h2>
<h3>Lat/Lon (Geographic) Coordinates should normally be avoided</h3>
<p>The Geographic Coordinate System is a spherical coordinate system, meaning that using them with Cartesian mathematics is technically incorrect. This demo uses them to increase the readability and is accurate enough to make the point. A better coordinate system would be Universal Transverse Mercator.</p>
<h3>Finding the corner coordinates</h3>
<p>One easy method to find the corner coordinates of an image is to use the command-line tool gdalinfo. For imagery which is ortho-rectified and contains the projection information, you can use the <a href="http://http://earthexplorer.usgs.gov">USGS EarthExplorer</a>. </p><div class="fragment"><div class="line">\f$&gt; gdalinfo N37W123.hgt</div><div class="line"></div><div class="line">   Driver: SRTMHGT/SRTMHGT File Format</div><div class="line">   Files: N37W123.hgt</div><div class="line">   Size is 3601, 3601</div><div class="line">   Coordinate System is:</div><div class="line">   GEOGCS[&quot;WGS 84&quot;,</div><div class="line">   DATUM[&quot;WGS_1984&quot;,</div><div class="line"></div><div class="line">   ... more output ...</div><div class="line"></div><div class="line">   Corner Coordinates:</div><div class="line">   Upper Left  (-123.0001389,  38.0001389) (123d 0&#39; 0.50&quot;W, 38d 0&#39; 0.50&quot;N)</div><div class="line">   Lower Left  (-123.0001389,  36.9998611) (123d 0&#39; 0.50&quot;W, 36d59&#39;59.50&quot;N)</div><div class="line">   Upper Right (-121.9998611,  38.0001389) (121d59&#39;59.50&quot;W, 38d 0&#39; 0.50&quot;N)</div><div class="line">   Lower Right (-121.9998611,  36.9998611) (121d59&#39;59.50&quot;W, 36d59&#39;59.50&quot;N)</div><div class="line">   Center      (-122.5000000,  37.5000000) (122d30&#39; 0.00&quot;W, 37d30&#39; 0.00&quot;N)</div><div class="line"></div><div class="line">    ... more output ...</div></div><!-- fragment --> <h2>Results </h2>
<p>Below is the output of the program. Use the first image as the input. For the DEM model, download the SRTM file located at the USGS here. <a href="http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N37W123.hgt.zip">http://dds.cr.usgs.gov/srtm/version2_1/SRTM1/Region_04/N37W123.hgt.zip</a></p>
<div class="image">
<img src="../../gdal_output.jpg" alt="gdal_output.jpg"/>
<div class="caption">
Input Image</div></div>
 <div class="image">
<img src="../../gdal_heat-map.jpg" alt="gdal_heat-map.jpg"/>
<div class="caption">
Heat Map</div></div>
 <div class="image">
<img src="../../gdal_flood-zone.jpg" alt="gdal_flood-zone.jpg"/>
<div class="caption">
Heat Map Overlay</div></div>
</div></div><!-- contents -->
<!-- HTML footer for doxygen 1.8.6-->
<!-- start footer part -->
<hr class="footer"/><address class="footer"><small>
Generated on Fri Apr 2 2021 11:36:34 for OpenCV by &#160;<a href="http://www.doxygen.org/index.html">
<img class="footer" src="../../doxygen.png" alt="doxygen"/>
</a> 1.8.13
</small></address>
<script type="text/javascript">
//<![CDATA[
addTutorialsButtons();
//]]>
</script>
</body>
</html>
